Symbiont ecology and its influence on the physiology and isotopic composition of a reef coral across space and time
CB Wall cbwall@hawaii.edu
9/19/2018
require("knitr")
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.align='center')
library(ggplot2)
library(plyr)
library(dplyr)
library(plotrix)
library(lubridate)
library(effects)
library(RCurl)
library(chron)
library(lubridate)
library(lme4)
library(emmeans)
library(multcomp)
library(devtools)
library(logihist)
library(popbio)
library(car)
library(lmerTest)
library(cowplot)
library(sjPlot)
library(sjmisc)Project background
Environmental
Sea level correction
Corals were collected across in summer and winter 2016 across several different days with each period. To correct for differences in the depth of coral colonies we used NOAA tide data to correct all depth measurements to Mean Sea Level (MSL). This approach was done in Innis et al. 2018 on a larger collection of the coral samples–which have been subset in the current study.
#########################
# Sea level correction
#########################
#### attach data files
data<-read.csv("data/mastersheet_PanKBAY.csv") # master file
qPCR.Innis<-read.csv("data/PanKBay_summer_qPCR.csv") # qPCR from summer (Innis et al. 2018)
JuneTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160801-20160831.csv")
DecemberTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20161201-20161222.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide, DecemberTide) # all MSL in meters
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
# use attr(as.POSIXlt(Tide$Time),"tzone") to confirm in PST
Tide<-Tide[, c(-5:-8)] #remove unnecessary columns
data$date.time<- as.POSIXct(paste(as.character(data$Date), as.character(data$Time.of.collection)),
format="%m/%d/%y %H:%M", tz="Pacific/Honolulu") # make sure time in HST for master
data$Time.r<-round_date(data$date.time,unit="6 minutes")
Tide$Time.r <- Tide$Time
data<-merge(data, Tide, by="Time.r", all.x=T)
data$newDepth <- data$Depth..m-data$TideHT # new depth in METERSLight parameters
DLI 2m all sites
This figure shows the daily integrated light (DLI mol photon m-2 d-1) at long-term light loggers at 2m depth near the locations where corals were collected in northern and southern sectors of Kāne’ohe Bay near shore (east) further off shore (west).
# DLI
##### all DLI for graphs
files <- list.files(path="data/environmental/temp and light/Jun_DecPAR/all DLI", pattern = "csv$", full.names = T)
tables <- lapply(files, read.csv, header = TRUE)
DLI<-do.call(rbind, tables)
DLI=DLI[-1] # remove junk column
DLI$Date<-as.Date(DLI$Date) # fix date
DLI<-DLI[order(DLI$Date),] # order by Date
# make a date sequence for entire study
all.date<-as.data.frame(seq(as.Date("2016-06-10"), as.Date("2017-01-12"), "days"))
colnames(all.date)="Date"
# merge the DLI data and the date sequence to make a complete df through time
DLI.3site<-merge(all.date, DLI, by="Date", all.x=T)
R10<-read.csv("data/environmental/temp and light/Jun_DecPAR/all DLI/Reef 10/DLI.Rf102016.csv")
R10<-R10[-1]; R10$Date<-as.Date(R10$Date)
DLI.4site<-merge(DLI.3site, R10, by="Date", all.x=T)
DLI.4site$month <- months(as.Date(DLI.4site$Date)) # makes a month column
DLI.4site<-DLI.4site[, c(1,6, 2:5)]
# write.csv(DLI.4site, "data/environmental/temp and light/Jun_DecPAR/All.DLI.csv")
##### Figure
All.DLI<-DLI.4site
reefcols=c("mediumseagreen", "dodgerblue", "salmon", "orchid")
par(mar=c(2,3.6,1,1.3), mgp=c(2,0.5,0))
layout(matrix(c(1,2,3,4), 2,2, byrow=TRUE))
plot(Rf44.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[1])
legend("topright", lty=1, col=reefcols[1], legend="Northwest", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
plot(Rf42.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[2])
legend("topright", lty=1, col=reefcols[2], legend="Northeast", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
plot(Rf10.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[3])
legend("topright", lty=1, col=reefcols[3], legend="Southwest", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
plot(HIMB.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[4])
legend("topright", lty=1, col=reefcols[4], legend="Southeast", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
Figure DLI at the four reef areas from summer to winter 2016 at 2m depths. Data at 2m was recorded from loggers
dev.copy(pdf, "figures/environmental/All.DLI.pdf", width=6.5, height=4)
dev.off()Model results testing differences in DLI at 2 m depth at all sites.
mod.light<-DLI.4site
Rf44DLI.mod<-mod.light[1:3]; Rf44DLI.mod$Location<-"Northwest";
colnames(Rf44DLI.mod)<-c("Date", "month", "DLI", "Location")
Rf42DLI.mod<-mod.light[,c(1:2,4)]; Rf42DLI.mod$Location<-"Northeast";
colnames(Rf42DLI.mod)<-c("Date", "month", "DLI", "Location")
HIMBDLI.mod<-mod.light[,c(1:2,5)]; HIMBDLI.mod$Location<-"Southeast";
colnames(HIMBDLI.mod)<-c("Date", "month", "DLI", "Location")
Rf10DLI.mod<-mod.light[,c(1:2,6)]; Rf10DLI.mod$Location<-"Southwest";
colnames(Rf10DLI.mod)<-c("Date", "month", "DLI", "Location")
All.DLI.mod<-rbind(Rf44DLI.mod, Rf42DLI.mod, HIMBDLI.mod, Rf10DLI.mod)
All.DLI.mod$Location<-as.factor(All.DLI.mod$Location)
All.DLI.mod$month<-as.factor(All.DLI.mod$month)
All.DLI.mod$Season<-ifelse(All.DLI.mod$month== "June" |
All.DLI.mod$month== "July" |
All.DLI.mod$month== "August" |
All.DLI.mod$month== "September", "dry season", "wet season")
All.DLI.mod$Season<-as.factor(All.DLI.mod$Season)
mod<-lmer(DLI~Location*Season + (1|Date), data=All.DLI.mod); anova(mod, type=2)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Location 4378.8 1459.58 3 529.66 134.674 < 2.2e-16 ***
## Season 1040.9 1040.91 1 210.48 96.043 < 2.2e-16 ***
## Location:Season 490.9 163.62 3 531.14 15.097 1.924e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(mod))
posthoc<-emmeans(mod, ~Location|Season)
CLD(posthoc, Letters=letters)## Season = dry season:
## Location emmean SE df lower.CL upper.CL .group
## Southwest 8.149681 0.4014748 504.55 7.360912 8.938449 a
## Northwest 10.829854 0.4450350 590.91 9.955811 11.703896 b
## Southeast 14.748648 0.4016675 503.81 13.959498 15.537797 c
## Northeast 14.997825 0.4601858 615.62 14.094100 15.901549 c
##
## Season = wet season:
## Location emmean SE df lower.CL upper.CL .group
## Southwest 3.901772 0.4473755 539.33 3.022959 4.780584 a
## Southeast 8.607241 0.4662620 551.22 7.691373 9.523109 b
## Northwest 9.128302 0.4662620 551.22 8.212435 10.044170 b
## Northeast 9.448888 0.4707121 559.55 8.524310 10.373467 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## wet season 7.771551 0.3348801 206.76 7.111333 8.431768 a
## dry season 12.181502 0.3071824 206.48 11.575885 12.787118 b
##
## Results are averaged over the levels of: Location
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05
Figure DLI model output
PAR and DLI at 3 depths
With Depth as a factor (<1m, 2m, 8m), we can produce a continuous graph of light (PAR or DLI) at 2m and use model intercepts to predict the light at the shallow and deep zones (ca. <1m and 8m, respectfully). Here, the model is \(lm(log.PAR\) ~ \(Depth.factor)\).
# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))
plot(y=All.PAR.df$Rf44.shall, x=All.PAR.df$timestamp, type="l", col="palegreen3",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="", main="Northwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.deep, x=All.PAR.df$timestamp, type="l", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
# Plot it!!
# Reef 42
plot(y=All.PAR.df$Rf42.shall, x=All.PAR.df$timestamp, type="l", col="mediumturquoise",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="", main="Northeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.deep, x=All.PAR.df$timestamp, type="l", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
# Plot it!!
# Reef 10
plot(y=All.PAR.df$Rf10.shall, x=All.PAR.df$timestamp, type="l", col="orangered",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.deep, x=All.PAR.df$timestamp, type="l", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
# Plot it!
# HIMB
plot(y=All.PAR.df$HIMB.shall, x=All.PAR.df$timestamp, type="l", col="mediumorchid2",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.deep, x=All.PAR.df$timestamp, type="l", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
Figure PAR at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations
dev.copy(pdf, "figures/environmental/PAR.all depths.pdf", height=7, width=8)
dev.off()#######
####### calculate DLI and compare
#######
df<-All.PAR.df
df$timestamp<-strptime(df$timestamp, format="%Y-%m-%d %H:%M:%S")
df$timestamp<-as.Date(df$timestamp, format="%Y-%m-%d") # change to DATE ONLY format to calulate range, means
df[is.na(df)] <- 0
df.split <- split(df, f=df$timestamp
< as.Date("2016-06-10", format="%Y-%m-%d")) # split df by date
df.dli<-aggregate(data.frame(Rf42.shall.DLI=df.split[[1]]$Rf42.shall*0.0864,
Rf42.mid.DLI=df.split[[1]]$Rf42.mid*0.0864,
Rf42.deep.DLI=df.split[[1]]$Rf42.deep*0.0864,
Rf44.shall.DLI=df.split[[1]]$Rf44.shall*0.0864,
Rf44.mid.DLI=df.split[[1]]$Rf44.mid*0.0864,
Rf44.deep.DLI=df.split[[1]]$Rf44.deep*0.0864,
HIMB.shall.DLI=df.split[[1]]$HIMB.shall*0.0864,
HIMB.mid.DLI=df.split[[1]]$HIMB.mid*0.0864,
HIMB.deep.DLI=df.split[[1]]$HIMB.deep*0.0864,
Rf10.shall.DLI=df.split[[1]]$Rf10.shall*0.0864,
Rf10.mid.DLI=df.split[[1]]$Rf10.mid*0.0864,
Rf10.deep.DLI=df.split[[1]]$Rf10.deep*0.0864),
by=list(Date=df.split[[1]]$timestamp), FUN=mean)
# make a date sequence for entire study
all.date.time<-as.data.frame(seq(
from=as.POSIXct("2016-06-10", tz="HST"),
to=as.POSIXct("2017-01-12", tz="HST"),
by="1 d"))
colnames(all.date.time)[1]<-"Date"
all.date.time$Date<-strptime(all.date.time$Date, format="%Y-%m-%d")
all.date.time$Date<-as.Date(all.date.time$Date, format="%Y-%m-%d")
df.dli<-merge(all.date.time, df.dli, by="Date", all.x=T)
df.dli[df.dli<=0] <- NA
#####################
# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))
plot(y=df.dli$Rf44.shall, x=df.dli$Date, type="h", col="palegreen3",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="", main="Northwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.deep, x=df.dli$Date, type="h", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
# Plot it!!
# Reef 42
plot(y=df.dli$Rf42.shall, x=df.dli$Date, type="h", col="mediumturquoise",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="", main="Northeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.deep, x=df.dli$Date, type="h", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
# Plot it!!
# Reef 10
plot(y=df.dli$Rf10.shall, x=df.dli$Date, type="h", col="orangered",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.deep, x=df.dli$Date, type="h", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
# Plot it!
# HIMB
plot(y=df.dli$HIMB.shall, x=df.dli$Date, type="h", col="mediumorchid2",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.deep, x=df.dli$Date, type="h", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
Figure DLI at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations
dev.copy(pdf, "figures/environmental/DLIcalc.all depths.pdf", height=7, width=8)
dev.off()DLI.summary<-read.csv("data/environmental/temp and light/Jun_DecPAR/DLI.3depth.summary.csv")
DLI.summary$Location<-factor(DLI.summary$Location, levels=c("NW", "NE", "SW", "SE"))
DLI.summary$Depth.zone<-factor(DLI.summary$Depth.zone, levels=c("shallow", "mid", "deep"))
#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
theme(text=element_text(size=10),
axis.line=element_blank(),
legend.justification=c(1,1), legend.position=c(1,1),
legend.background = element_rect("transparent", colour=NA),
legend.text=element_text(size=10),
legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=0.8),
aspect.ratio=1,
axis.ticks = element_line(colour = 'black', size = 0.4),
axis.ticks.length=unit(0.25, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
legend.key.size = unit(0.4, "cm"))
pd <- position_dodge(0.7) #offset for error bars and columns
#####
Location.label<-c("NW", "NE", "SW", "SE")
Depth.label<-c("<1m", "2m", "8m")
plot_colors2<-c("lightskyblue", "dodgerblue", "deepskyblue4")
plot_stat<-c("slategray", "dodgerblue2", "steelblue1", "slategray", "dodgerblue2", "steelblue1",
"slategray", "dodgerblue2", "steelblue1","slategray", "dodgerblue2", "steelblue1")
# DLI by site and depth
Fig.DLI<-ggplot(DLI.summary, aes(x=Location, y=DLI.mean, fill=Depth.zone)) +
geom_bar(stat="identity", col=plot_stat, position = pd, width=0.7, size=0.3) +
geom_errorbar(aes(ymin=DLI.mean-DLI.se, ymax=DLI.mean+DLI.se),
size=0.4, width=0, position=pd, col=plot_stat) +
xlab("Locations") +
ylab(expression(paste("DLI", ~(mol~photons~m^-2~d^-1), sep=""))) +
scale_x_discrete(labels=Location.label) +
scale_y_continuous(expand=c(0,0), limits=c(0, 35)) +
scale_fill_manual(values=alpha(c(plot_colors2), 0.6),
labels=Depth.label)+ Fig.formatting; Fig.DLI
Figure Daily Light Integral (DLI) for four reef locations at three depths. DLI estimated from light attenuation coefficents using depth as a fixed effect against baseline loggers at 2 m depth
ggsave(file="figures/environmental/DLI.bar.pdf", height=4, width=4)Attenuation and Light-at-Depth
A periodic deployment of loggers (October and November) at ca. <1m and ca. 7m were used to model the relationship between 2m light and light at other depths.
Using the logger data, the difference in light at 1 and 7m was calculated relative to 2m (i.e, PAR baseline (2m) - PAR other), then the difference in depth at 1 and 7m was calculated realtive to logger depth i.e, depth logger other - depth logger (2m baseline). A no-intercept model was used to fit difference in log(DLI) and depth (i.e, log(delta-DLI) and delta-depth) with depth as a continous/numeric variable to determine kd: \(lm(delta.log.DLI\)~\(delta.Depth + 0)\).
Site-specific kd values were used to generate a continuous variable of “light at depth”. At each time period (summer or winter) the daily light integral was averaged over 3 months: summer 2016 (June, July, August) and winter 2016 (November, December, January). Discrete monthly seasonal-mean DLI means was used to generate an estimate of light at depth for corals at each site, as this integrates the site and seasoanl variance.
Light at Depth was calculated following Beer-Lambert: \(Ez{_d} = Ez_{_0}e^{(-k_d *d)}\), where \(d\) is depth, \(k_d\) is the attenuation coefficient for each site, \(Ez{_d}\) is light at depth \(d\) and \(Ez{_0}\) is measured light. Calculations were relative to light at 2m depth (i.e., \(Ez_{_0}\) was the light at 2m depth) and \(d\) was the difference between coral depth relative to depth of the logger (i.e., 2m). light reaching depth d and Ez(0) is the incident light at the surface.
########################## # LIGHT AT DEPTH
#############
######
# using new depth and the site-specific kd (light attenuation) determine approximate "light at depth" for each colony sample.
######
### attach necessary files
logger.depths<-read.csv("data/environmental/temp and light/light.logger.depths.csv") # depths for loggers
kd.all<-read.csv("data/environmental/kd.all.csv") # kds for each site
# Seasonal DLI used for "period of collection" light levels
month.DLI<-read.csv("data/environmental/temp and light/Jun_DecPAR/All.DLI_long.csv")
# corals collected in June-July-August use summer time DLI for these months as indicator of average DLI
summer.DLI<-month.DLI[(month.DLI$Month=="June" | month.DLI$Month=="July" | month.DLI$Month=="August"),]
winter.DLI<-month.DLI[(month.DLI$Month=="November" | month.DLI$Month=="December" |month.DLI$Month=="January"),]
# summer mean and SE dataframe
sum.mean<-aggregate(DLI~Site, summer.DLI, mean)
sum.SE<-aggregate(DLI~Site, summer.DLI, std.error)
sum.light.df<-cbind(sum.mean, sum.SE[2])
colnames(sum.light.df)=c("Site", "mean.DLI", "se.DLI")
sum.light.df$Season<-as.factor("summer")
# winter mean and SE dataframe
wint.mean<-aggregate(DLI~Site, winter.DLI, mean)
wint.SE<-aggregate(DLI~Site, winter.DLI, std.error)
wint.light.df<-cbind(wint.mean, wint.SE[2])
colnames(wint.light.df)=c("Site", "mean.DLI", "se.DLI")
wint.light.df$Season<-as.factor("winter")
season.DLI<-rbind(sum.light.df[,c(4,1,2:3)], wint.light.df[,c(4,1,2:3)]) # compiled means for DLI at ~2m
season.DLI$Location<- as.factor(c("SE", "SW", "NE", "NW", "SE", "SW", "NE", "NW")); season.DLI<-season.DLI[,c(1,5,2:4)]
write.csv(season.DLI, "data/environmental/season.DLI.csv")
#######################################
### make new dataframe for calculations
df.light<- data[, c("Season", "Reef", "Location", "newDepth")]
# make a column of "depth differences" relative to where ~2m logger was deployed
df.light$depth.diff<-ifelse(df.light$Reef=="F1-R46", df.light$newDepth-1.8288,
ifelse(df.light$Reef=="F8-R10", df.light$newDepth-1.8288,
ifelse(df.light$Reef=="HIMB", df.light$newDepth-1.8288,
df.light$newDepth-2.4384))) # last statement for remaining site (Rf42)
# make a column for sample-specific light at depth (estimate) based on kd
# follow: #with 2m as baseline# E(depth) = E(2m)*exp(-k_d*(depth - 2m))
# so that: light at depth = DLI at mid.depth * exp (-kd *(delta shallow-deep in m))
df.light$Light<-ifelse(df.light$Reef=="HIMB" & df.light$Season=="summer",
season.DLI$mean.DLI[1]*exp(-kd.all$kd[1]*df.light$depth.diff), # summer HIMB
ifelse(df.light$Reef=="F8-R10" & df.light$Season=="summer",
season.DLI$mean.DLI[2]*exp(-kd.all$kd[2]*df.light$depth.diff), # summer R10
ifelse(df.light$Reef=="R42" & df.light$Season=="summer",
season.DLI$mean.DLI[3]*exp(-kd.all$kd[3]*df.light$depth.diff), # summer R42
ifelse(df.light$Reef=="F1-R46" & df.light$Season=="summer",
season.DLI$mean.DLI[4]*exp(-kd.all$kd[4]*df.light$depth.diff), # summer R46
ifelse(df.light$Reef=="HIMB" & df.light$Season=="winter",
season.DLI$mean.DLI[5]*exp(-kd.all$kd[1]*df.light$depth.diff), # winter HIMB
ifelse(df.light$Reef=="F8-R10" & df.light$Season=="winter",
season.DLI$mean.DLI[6]*exp(-kd.all$kd[2]*df.light$depth.diff), # winter R10
ifelse(df.light$Reef=="R42" & df.light$Season=="winter",
season.DLI$mean.DLI[7]*exp(-kd.all$kd[3]*df.light$depth.diff), # winter R42
season.DLI$mean.DLI[8]*exp(-kd.all$kd[4]*df.light$depth.diff)) # winter R46
))))))
###### plot of light x depth by season
df.light$Reef <- factor(df.light$Reef, levels = c("F1-R46", "R42", "F8-R10", "HIMB"))
df.light$Location <- factor(df.light$Location, levels=c("NW", "NE", "SW", "SE"))
plot.by.sites=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Location), alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons")
plot.by.site=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Season), alpha=0.3, position = 'stack') +
scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") + facet_wrap(~Location, scales="free") + scale_fill_manual(values=c("darkorange1", "dodgerblue1"))
######
# can loop as a list
p=ggplot(df.light, aes(Light, fill=Season)) + geom_density(alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons")
plots=dlply(df.light, .(Location), function(x) p %+% x + facet_wrap(~Location))
###### plot of light x depth by season with exponential curve fitting
Sum<-df.light[(df.light$Season=="summer"),]
Win<-df.light[(df.light$Season=="winter"),]
plot(Light~newDepth, Sum, col="coral", pch=16, xlab="Depth (m)", ylab="DLI at coral",
ylim=c(0, 30), xlim=c(0, 10))
summod<-lm(log(Light)~newDepth, Sum)
curve(exp(coef(summod)["(Intercept)"]+coef(summod)["newDepth"]*x), add=TRUE, col="coral", lwd=2)
par(new=T)
plot(Light~newDepth, Win, col="lightblue2", pch=16, xaxt="n", yaxt="n",
xlab="", ylab="", ylim=c(0, 30), xlim=c(0, 10))
wintmod<-lm(log(Light)~newDepth, Win)
curve(exp(coef(wintmod)["(Intercept)"]+coef(wintmod)["newDepth"]*x), add=TRUE, col="lightblue2", lwd=2)
legend("topright", legend=c("summer", "winter"), col=c("coral", "lightblue2"), lty=1, lwd=1, pch=16, bty="n")
Figure Relationship between mean seasonal light availability (mean daily light integral) at depth of coral fragment collection
Dissolved nutrients
Dissolved inorganic nutrients in seawater were quantified on water samples (ca. 100 ml) taken from each site and filtered through a GF/F filter and stored in acid-washed bottles and frozen at -20 °C until analyzes. Dissolved inorganic nutrients—ammonium (NH4+), nitrate + nitrite (NO3- + NO2-), phosphate (PO43-), and silicate (Si(OH)4)—were measured at University of Hawai‘i at Mānoa SOEST Laboratory for Analytical Biogeochemistry using a Seal Analytical AA3 HR nutrient autoanalyzer and expressed as μmol l-1.
Silicate showed no effects (p > 0.323). Phosphate increased in the winter (p=0.046) most notably at northern sites. N+N was higher in northern sites relative to southern sites and increased during the winter at all sites compared to the summer. Ammonium increased in the winter as well (p <0.001).
#########################
#########################
# nutrients
nutr<-read.csv("data/environmental/PanKBay_nutrients.csv")
#models
mod<-lm(silicate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(phosphate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(N.N~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(ammonium~Location+Date, data=nutr); Anova(mod, type=2)
nutr$Date<-mdy(nutr$Date) # corrects date format
dates<-cbind("10-Aug '16", "20-Dec '16")
RfHIMB<-nutr[(nutr$Reef=="HIMB"), ]
Rf42<-nutr[(nutr$Reef=="R42"), ]
Rf46<-nutr[(nutr$Reef=="F1-46"), ]
RF10<-nutr[(nutr$Reef=="F8-R10"), ]
Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Locations<-c("NW", "NE", "SW", "SE")
plot_colors<-c("mediumseagreen", "dodgerblue", "salmon", "orchid")
###############
# Phosphate
# par(mfrow=c(1,1), mar=c(2,4,1,1), mgp=c(2,0.5,0))
# figure layout
layout(matrix(c(1,2,3,4), nrow=1, byrow=TRUE))
par(mar=c(5,5,7,1))
###############
# Silicate
plot(y=Rf46$silicate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("silicate"~(mu*mol~L^-1), sep="")), ylim=c(0, 15), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46 / NW
#plots Reef 42/ NE
with(Rf46, lines(Rf42$silicate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10/ SW
with(Rf46, lines(RF10$silicate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB/ SE
with(RfHIMB, lines(RfHIMB$silicate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
###############
# Phosphate
plot(y=Rf46$phosphate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("phosphate"~(mu*mol~L^-1), sep="")), ylim=c(0, 0.25), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Rf46 /NW
#plots Reef 42 / NE
with(Rf46, lines(Rf42$phosphate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10 /SW
with(Rf46, lines(RF10$phosphate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB / SE
with(RfHIMB, lines(RfHIMB$phosphate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
###############
# N+N
plot(y=Rf46$N.N, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("nitrate+nitrite"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46/ NW
#plots Reef 42/ NE
with(Rf46, lines(Rf42$N.N, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10 / SE
with(Rf46, lines(RF10$N.N, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB / SW
with(RfHIMB, lines(RfHIMB$N.N, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
###############
# Ammonium
plot(y=Rf46$ammonium, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("ammonium"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46/ NW
#plots Reef 42 /NE
with(Rf46, lines(Rf42$ammonium, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10/ SW
with(Rf46, lines(RF10$ammonium, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB /SE
with(RfHIMB, lines(RfHIMB$ammonium, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
legend("topleft", legend=Locations, col=plot_colors[1:4], lwd=1, pch=19, lty=2, cex=0.8, bty="n", x.intersp=0.2, y.intersp=1.2, inset=c(0.01, 0), seg.len=1.5)
dev.copy(pdf, "figures/environmental/all.nutrients.pdf", encod="MacRoman", height=3, width=7)
dev.off()Suspended particulates
Gravimetric SPM
Total (organic + inorganic) suspended particulate matter (SPM) was quantified by collecting and filtering 0.5 L of seawater from each location during each season. Seawater was sieved at 243 μm to remove large debris and filtered onto a pre-combusted (4h, 450°C) GF/F filter (0.7 μm). Salts were removed by light rinsing with ddH2O. Filters were dried at 60°C for 24h and re-weighed to closest 0.001 g to obtain total SPM. Masses were normalized to a liter of water.
###### SPM
spm<-read.csv("data/environmental/PanKBay_SPM.csv")
spm$SPM..g.l<-spm$SPM..g.l*1000 # change to mg
spm$Date<-factor(spm$Date, levels=c("8/10/16", "12/20/16"))
colnames(spm)[4]="SPM.mg"
spm$Location<-ordered(spm$Location, c("NW", "NE", "SW", "SE"))
mod<-lm(SPM.mg~Location+Date, data=spm); Anova(mod, type=2)## Anova Table (Type II tests)
##
## Response: SPM.mg
## Sum Sq Df F value Pr(>F)
## Location 8.8382 3 10.3062 0.0003025 ***
## Date 1.4114 1 4.9373 0.0386252 *
## Residuals 5.4312 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod, ~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## NE 0.7966668 0.2182711 19 0.3398203 1.253513 a
## NW 1.6450002 0.2182711 19 1.1881536 2.101847 ab
## SE 2.2666668 0.2182711 19 1.8098203 2.723513 b
## SW 2.2850002 0.2182711 19 1.8281536 2.741847 b
##
## Results are averaged over the levels of: Date
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05plot(allEffects(mod), par.strip.text=list(cex=0.7))
rm(mean) # this will remove 'mean' from library and let the function below proceed.
mean<-aggregate(SPM.mg~Date + Location+Reef, spm, mean)
se<-aggregate(SPM.mg~Date+Location+Reef, spm, std.error)
spm.df<-cbind(mean, se[c(0,4)])
colnames(spm.df)=c("Date", "Location", "Reef", "SPM", "se")
spm.df$Date<-ordered(spm.df$Date, c("8/10/16", "12/20/16"))
spm.df$Reef<-ordered(spm.df$Reef, c("F1-46", "R42", "F8-10", "HIMB"))
spm.df$Location<-ordered(spm.df$Location, c("NW", "NE", "SW", "SE"))
Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Locations<-c("NW", "NE", "SW", "SE")
Date<-c("10-Aug '16", "20-Dec '16")
plot_colors2<-c("mediumturquoise", "skyblue3", "lightcoral", "plum")
#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
theme(text=element_text(size=10),
axis.line=element_blank(),
legend.justification=c(1,1), legend.position=c(1,1),
legend.background = element_rect("transparent", colour=NA),
legend.text=element_text(size=10),
legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=0.8),
aspect.ratio=1,
axis.ticks = element_line(colour = 'black', size = 0.4),
axis.ticks.length=unit(0.25, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
legend.key.size = unit(0.4, "cm"))
pd <- position_dodge(0.7) #offset for error bars and columns
#####
# SPM (mg/l)
Fig.SPM<-ggplot(spm.df, aes(x=Date, y=SPM, fill=Location)) +
geom_bar(colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3", "mediumturquoise"),
stat="identity", position = pd, width=0.7, size=0.3) +
geom_errorbar(aes(ymin=SPM-se, ymax=SPM+se),
colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3","mediumturquoise"),
size=0.4, width=0, position=pd) +
xlab("Sampling points") +
ylab(expression(paste("SPM", ~(mg~L^-1), sep=""))) +
scale_x_discrete(labels=Date) +
scale_y_continuous(expand=c(0,0), limits=c(0, 5)) +
scale_fill_manual(values=alpha(c(plot_colors2), 0.5),
labels=Locations)+ Fig.formatting; Fig.SPM
ggsave(file="figures/environmental/SPM.pdf", height=4, width=4)Isotopes of SPM/plankton
Plankton tows and seawater collections to parameterized particulates and plankters as potential heterotrophic end members available to reef corals. Sampling was done at 4 sites where corals were collected [North: (Reef 42, fringe-Reef 46), and South: (HIMB, fringe-Reef 10)], as well as two reefs in central region of the bay where corals were not collected (Central: Reef 25, fringe-Reef 25) to increase spatial resolution of suspended particulate isotope sample sizes. Using size-fractioned materials collected in seawater and plankton tows, we examined the spatial (among reef sites) and temporal patterns (among seasons) in stable isotope values of size-fractioned end members. Carbon and nitrogen isotope values among the 6 reef sites were not significantly different (p > 0.140), season had marginal effects (p > 0.049), whereas fraction influenced isotope values significantly (p< 0.001). Therefore, isotope values were pooled among reefs and seasons to best interpret size-fraction isotope values. Samples were not acidified prior to analysis as to not affect nitrogen isotope values.
SWiso<-read.csv("data/isotopes_SW_all times.csv")
mod<-(lm(d13C~Location+Season+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no location effect## Anova Table (Type II tests)
##
## Response: d13C
## Sum Sq Df F value Pr(>F)
## Location 17.926 5 1.3418 0.2626013
## Season 7.921 1 2.9645 0.0914188 .
## SW.fraction..um 76.419 4 7.1504 0.0001286 ***
## Residuals 130.920 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod<-(lm(d15N~Location+Season+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no location effect## Anova Table (Type II tests)
##
## Response: d15N
## Sum Sq Df F value Pr(>F)
## Location 2.3255 5 1.7291 0.14557
## Season 1.0935 1 4.0653 0.04927 *
## SW.fraction..um 30.3773 4 28.2335 3.454e-12 ***
## Residuals 13.1802 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## summer 6.25 0.09468949 49 6.059714 6.440286 a
## winter 6.52 0.09468949 49 6.329714 6.710286 b
##
## Results are averaged over the levels of: Location, SW.fraction..um
## Confidence level used: 0.95
## significance level used: alpha = 0.05SWiso$Reef<-factor(SWiso$Reef, levels=c("F1-46", "R42", "F2-R25", "R25", "F8-R10", "HIMB"))
SWiso$Location<-factor(SWiso$Location, levels=c("NW", "NE", "CW", "CE", "SW", "SE"))
SWiso$SW.fraction..um<-factor(SWiso$SW.fraction..um, levels=c(">243", "<243", "100-243", "10-100", "0-10"))
winter.data<-SWiso[(SWiso$Season=="winter"),]
summer.data<-SWiso[(SWiso$Season=="summer"),]These box plots show the seasonal isotope values for carbon and nitrogen according to size fractions.
# plots of SW isotopes by Season-- first showing by size, then by site
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))
######### Sizes
# Summer size d13C
plot(d13C~SW.fraction..um, data=summer.data, ylim=c(-25,-10),
main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
# Winter size d13C
plot(d13C~SW.fraction..um, data=winter.data, ylim=c(-25,-10),
main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
# Summer size d15N
plot(d15N~SW.fraction..um, data=summer.data, ylim=c(3,10),
main=expression(paste("summer"~ delta^{15}, "N")), col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
xlab=expression(paste("Size Fraction (", mu,m,")")),
ylab=expression(paste(delta^{15}, "N (‰, air)")))
# Winter size d15N
plot(d15N~SW.fraction..um, data=winter.data, ylim=c(3,10),
main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
Figure Size-fractioned sample isotope values stratified by seasons (summer and winter)
These figures show the same size-fractioned isotope data as above, but stratified by location and pooled across seasons.
######## Sites
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))
# Summer site d13C
plot(d13C~Location, data=summer.data, ylim=c(-25,-10),
main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")
# Winter site d13C
plot(d13C~Location, data=winter.data, ylim=c(-25,-10),
main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")
# Summer site d15N
plot(d15N~Location, data=summer.data, ylim=c(3,10),
main=expression(paste("summer"~ delta^{15}, "N")), col="salmon", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab="Reef Sites")
# Winter site d15N
plot(d15N~Location, data=winter.data, ylim=c(3,10),
main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
col="salmon", cex.axis=0.7, cex.main=1, xlab="Reef Sites")
Figure Isotope values in seawater particles stratified by sites (northwest to southeast)
Pooling the isotope data across sites and seasons (where effects were absent or negligible), this planktonic foodweb for carbon and nitrogen illustrates differences in the food sources and their isotopic composition.
mod<-lm(d13C~SW.fraction..um, data=SWiso); Anova(mod, type=2)
posthoc<-emmeans(mod, ~SW.fraction..um)
CLD(posthoc, Letters=letters)
mod<-lm(d15N~SW.fraction..um, data=SWiso); Anova(mod, type=2)
posthoc<-emmeans(mod, ~SW.fraction..um)
CLD(posthoc, Letters=letters)
rm(mean)
####
#### making scatter for d15N and d13C, pooled across seasons and sites
mix.N.mean<-aggregate(d15N~SW.fraction..um, data=SWiso, mean)
mix.N.se<-aggregate(d15N~SW.fraction..um, data=SWiso, std.error)
mix.C.mean<-aggregate(d13C~SW.fraction..um, data=SWiso, mean)
mix.C.se<-aggregate(d13C~SW.fraction..um, data=SWiso, std.error)
mix.data<-cbind(mix.N.mean, mix.C.mean[c(2,0)], mix.N.se[c(2,0)], mix.C.se[c(2,0)]); colnames(mix.data)=c("fraction", "d15N", "d13C", "d15N.se", "d13C.se")
colors=c("#FF6A6A", "#00B2EE", "#FFB90F", "#3CB371", "#8B7500")
op<-par(mfrow = c(1,1), mar=c(5,4,1,5),xpd=TRUE, pty="sq")
size.labels=c(expression(paste("> 243"~mu,m)), expression(paste("< 243"~mu,m)), expression(paste("100-243"~mu,m)), expression(paste("10-100"~mu,m)), expression(paste("< 10"~mu,m)))
#### make the plot
plot(d15N~d13C, data=mix.data, type="n", ylim=c(4.5,8), xlim=c(-23,-17), tck=-0.03, cex.axis=0.7, cex.lab=0.7,
ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab=expression(paste(delta^{13}, "C (‰, V-PDB)")))
legend("topleft", inset=c(0.05,0.0), legend=size.labels, col=colors, pch=19, cex=0.6, bty="n", x.intersp=1.8, y.intersp = 1.4)
arrows(mix.data$d13C-mix.data$d13C.se, mix.data$d15N, mix.data$d13C+mix.data$d13C.se, mix.data$d15N, length=0.0, angle=90, code=3, lwd=0.7, col=colors)
arrows(mix.data$d13C, mix.data$d15N-mix.data$d15N.se, mix.data$d13C, mix.data$d15N+mix.data$d15N.se, length=0.0, angle=90, code=3, lwd=0.7, col=colors)
points(d15N~d13C, data=mix.data, pch=19, cex=0.8, col=colors)
Figure Isotope values for size-fractioned seawater particles and plankters pooled across seasons and reef sites
dev.copy(pdf, "figures/environmental/iso.sources.KBay.pdf", encod="MacRoman", height=4, width=4)
dev.off()Biology
################################
### Biological responses
################################
#data : this is the master file
# add in light at depth column from df.light dataframe
data$Light<-df.light$Light
##### produce a categorical depth bin ####
depth<-data$newDepth
data$depth.bin<-factor(ifelse(depth<2, "<2m", ifelse(depth >2 & depth <4, "2-4m", ifelse(depth >4 & depth <6, "4-6m", ">6m"))), levels=c("<2m", "2-4m", "4-6m", ">6m"))
aggregate(Sample.ID~depth.bin+Season+Location, data, length)
data$depth.bin.small<-factor(ifelse(depth<4, "<4m", ">4m"), levels= c("<4m", ">4m"))
################################################
# calculate, normalized dependent variables
################################################
str(data)
data$cells.ml<-as.numeric(data$cells.ml)
# helpful shorthand
SA<-data$surface.area # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml
# AFDW.mg. == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$biomass<- (data$mg.biomass.ml*blastate)/SA
# Symbiodinium.cells. == cell.ml * blastate / SA
data$zoox<- (data$cells.ml*blastate)/SA
# total chlorophyll == ug.chl.a.ml * blastate + ug.chl.c2.ml * blastate / SA
data$chltot<-(data$ug.chl.a.ml)+(data$ug.chl.c2.ml)*blastate/SA
# pg.chlorophyll.a..cell + pg.chlorophyll.c2..cell == ug.chltot.ml * 10^6 / cells.ml
data$chlcell<- (data$ug.chl.a.ml*10^6+data$ug.chl.c2.ml*10^6)/data$cells.mlqPCR
####################
# qPCR
#########
# qPCR
# Use steponeR to import data and calculate proporation of C and D symbionts
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path="data/qPCR", pattern = "txt$", full.names = T); Mcap.plates
Mcap <- steponeR(files=Mcap.plates, delim="\t",
target.ratios=c("C.D"),
fluor.norm=list(C=2.26827, D=0),
copy.number=list(C=33, D=3),
ploidy=list(C=1, D=1),
extract=list(C=0.813, D=0.813))
Mcap <- Mcap$result
head(Mcap)
# remove +/-control
Mcap <- Mcap[grep("+C52", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("H2O", Mcap$Sample.Name, fixed=T, invert = T), ]
# to remove any early-amplification CT noise
Mcap$C.CT.mean[which(Mcap$C.CT.mean < 15)] <- 0
#Remove failed samples, i.e., those where either C or D were NOT found in both reps
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
# replace CT means with 'NA' as zero
Mcap$C.CT.mean[is.na(Mcap$C.CT.mean)] <-0
Mcap$D.CT.mean[is.na(Mcap$D.CT.mean)] <-0
Mcap$C.D[is.na(Mcap$C.D)] <- 1 # sets all infinity (= 100% C) to 1.0
# caluclate proportion C and proprtion D where C and D are both present
Mcap$propC<- Mcap$C.D / (Mcap$C.D + 1)
Mcap$propD<- 1 / (Mcap$C.D + 1)
# where C and D are not cooccuring...
# if C.D = 1 = 100% C, make 'PropC' = 1 and 'PropD' = 0
# if C.D = 0 = 100% D, make 'PropD' = 1 and 'PropC' = 0
Mcap$propC[which(Mcap$C.D==1)] <- 1
Mcap$propD[which(Mcap$propC==1)] <- 0
Mcap$propD[which(Mcap$C.D==0)] <- 1
# calculate FOUR COMMUNITY categories: C, C>D, D>C, D
Mcap$Mix <- factor(ifelse(Mcap$propC > Mcap$propD, ifelse(Mcap$propD!= 0, "CD", "C"), ifelse(Mcap$propD > Mcap$propC, ifelse(Mcap$propC!=0, "DC", "D"), NA)), levels=c("C", "CD", "DC", "D"))
# Identify SINGLE dominant symbiont clade: C or D
Mcap$Dom <- factor(substr(as.character(Mcap$Mix), 1, 1))
# Set zeros to NA to facilitate log transformation
Mcap$propC[which(Mcap$propC==0)] <- NA
Mcap$propD[which(Mcap$propD==0)] <- NA
######## look for duplicates in dataset by year and type of event (bleach/recover)
Mcap[duplicated(Mcap$Sample.Name), ] ## duplicates
# remove duplicated
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_05" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_06" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_01" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_02" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_03" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_13" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_14" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate3.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_11" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
# parse Sample ID and Site from "Site.Name"
Mcap<-cbind(Mcap, colsplit(Mcap$Sample.Name, pattern= "_", c("Reef", "Sample.ID")))
Mcap$Season<-as.factor("winter")
Mcap$Reef<-as.factor(Mcap$Reef)
Mcap$Reef<-revalue(Mcap$Reef, c("R10"="F8-R10", "R46"="F1-R46")) # rename factor levels
# make new factors for bay region and reef type
Mcap$Bay.region <- ifelse(Mcap$Reef=="R42" | Mcap$Reef=="F1-R46", "northern", "southern")
Mcap$Reef.type <- ifelse(Mcap$Reef=="R42" | Mcap$Reef=="HIMB", "patch", "fringe")
Mcap$Location <- ifelse(Mcap$Reef=="F1-R46", "NW", ifelse(Mcap$Reef=="R42", "NE",
ifelse(Mcap$Reef=="F8-R10", "SW", "SE")))
### reorder columns and finish
Mcap<-Mcap[, c(17,15,20,19,18,16,1:14)] # reordered to match masterdata, and finish
### structure winter and summer qPCR dataframes to have same columns and combine dataframe
qPCR.winter<-Mcap[ , (names(Mcap) %in%
c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))]
qPCR.summer<-qPCR.Innis[ , (names(qPCR.Innis) %in%
c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))]
# merge qPCR files
qPCR.all<-rbind(qPCR.winter, qPCR.summer)
# add to master data
data.all<-merge(data, qPCR.all, by=c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID"), all.x=T)
###### remove columns no longer needed, update "Depth" to be tide-corrected depth )= newDepth
data.trim<-data.all[ , !(names(data.all) %in% c("total.blastate.ml", "Date", "Time.of.collection", "Depth..m", "Time.r", "surface.area.cm2", "cells.ml", "ug.chl.a.ml", "ug.chl.c2.ml", "mg.biomass.ml", "host..mass.mg", "host..ugN", "host..ugC", "symb..mass.mg", "symb..ugN", "symb..ugC", "stationId", "datum", "TimeUTC", "TideHT", "Time"))]
data.trim$symb..C.N[data.trim$symb..C.N>=12.520270]=NA # set this outlier to NAqPCR figure The distribution of C and D symbiont types is modeled using a logistic regression with dominant symbiont community (C vs. D) as a response and Depth (m), Season, and reef Location as a predictor. A generalized linear model with a binomial distribution and logit link function is used. Model selection (by model AIC values) between a fully crossed model (crossed) and additive model (add) were compared. The additive mode best fit the data; main effects were tested in anova with a Chi-square test.
## qPCR figures
#####
model.data<-data.trim[,c(17, 1:5, 18:25, 7:16, 26:29)]
#########################
#########################
#########################
# Plots with DEPTH as the predictor
# qPCR and symbionts
#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Reef <- factor(Symb$Reef, levels = c("F1-R46", "R42", "F8-R10", "HIMB")) # set levels
Symb$Location <- factor(Symb$Location, levels = c("NW", "NE", "SW", "SE")) # set levels
Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add <-glm(Dominant~newDepth+Season+Location, family = "binomial", data = Symb)
crossed <-glm(Dominant~newDepth*Season+Location, family = "binomial", data = Symb)
AIC(add, crossed) # additive model best## df AIC
## add 6 122.1298
## crossed 7 119.7058anova(crossed, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 118 140.51
## newDepth 1 25.9322 117 114.58 3.536e-07 ***
## Season 1 0.1044 116 114.48 0.74660
## Location 3 4.3477 113 110.13 0.22628
## newDepth:Season 1 4.4240 112 105.71 0.03544 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results.all<-crossed
#summary(results.all)
plot(allEffects(results.all), par.strip.text=list(cex=0.6))
results.all.mod <-glm(Dominant~newDepth, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum=glm(Dominant~newDepth, family = "binomial", data = sum.dat)
anova(results.sum, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 59 67.48
## newDepth 1 21.71 58 45.77 3.171e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))
###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win=glm(Dominant~newDepth, family = "binomial", data = wint.dat)
anova(results.win, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 58 72.583
## newDepth 1 7.1278 57 65.455 0.00759 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#summary(results.win)
fitted.win <- predict(results.win, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))Logistic regression of the probability of Montipora capitata have a clade D-dominant symbiont community (probability of 1.0) versus one dominated by clade C symbiont (probability of 0). Both seasons show similar patterns, although the probability of corals being clade D-dominant increases in winter months as a function of Depth relative to summer months.
######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$newDepth, xlab="Depth (m)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,10), ylim=c(0.0, 1.0), cex=0.5, cex.lab=0.8)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,10), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5, cex.lab=0.8)
lines(fitted.all ~ seq(0,10,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,10,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,10,0.1), col="lightskyblue", lwd=1)
legend("topright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.8, y.intersp=2, cex=0.5, inset=c(-0.005, -0.01), seg.len=4)
Figure Probability of Montipora capitata hosting clade D symbionts across a depth gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data
dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_depth.pdf", encod="MacRoman", height=4, width=4)
dev.off()Seen another way–Histograms. By separating these data into 2 histograms the probability of clade D- (1.0) and clade C-dominance (0.0) can be visualized. Here it is clear that there is a greater probability of finding corals with clade D dominant communites at shallow depths, whereas in winter there is a slightly greater probability of finding clade D in corals with increasing depth.
#######
#######
par(mfrow=c(1,2))
Dom <- subset(Symb, !is.na(newDepth) & !is.na(Dominant))
Dom.sum<-Dom[(Dom$Season=="summer"),]
Dom.win<-Dom[(Dom$Season=="winter"),]
logi.hist.plot(Dom.sum$newDepth, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)
logi.hist.plot(Dom.win$newDepth, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C D", line = 0.5, cex = 0.8)
Figure Histograms of probability of Montipora capitata hosting clade D symbionts across a depth gradient among seasons. Lines represent model fit to data by logistic regression
dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()#
#########################
#########################
# this code runs the same as above, except using colony LIGHT-AT-DEPTH as a predictor instead of depth
# Depth as a predictor
# qPCR and symbionts
#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add<-glm(Dominant~Light+Season+Location, family = "binomial", data = Symb)
crossed<-glm(Dominant~Light*Season*Location, family = "binomial", data = Symb)
#AIC(add, crossed)
anova(add, test = "Chisq")
#summary(add)
results.all<-add
results.all.mod<-glm(Dominant~Light, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(Light=seq(0,25,0.1)), type = "response")
###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum<-glm(Dominant~Light, family = "binomial", data = sum.dat)
#anova(results.sum, test = "Chisq")
#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))
###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win<-glm(Dominant~Light, family = "binomial", data = wint.dat)
#anova(results.win, test = "Chisq")
#summary(results.win)
fitted.win <- predict(results.win, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))
######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$Light, xlab="Light (DLI)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,25), ylim=c(0.0, 1.0), cex=0.5, cex.lab=0.8)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,25), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5)
lines(fitted.all ~ seq(0,25,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,25,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,25,0.1), col="lightskyblue", lwd=1)
legend("bottomright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.7, y.intersp=2, cex=0.5, inset=c(0.1, 0.15), seg.len=4)
dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_light.pdf", encod="MacRoman", height=4, width=4)
dev.off()
#######
#######
par(mfrow=c(1,2))
Dom <- subset(Symb, !is.na(Light) & !is.na(Dominant))
Dom.sum<-Dom[(Dom$Season=="summer"),]
Dom.win<-Dom[(Dom$Season=="winter"),]
logi.hist.plot(Dom.sum$Light, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Light-at-depth (DLI)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)
logi.hist.plot(Dom.win$Light, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Light-at-depth (DLI)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C D", line = 0.5, cex = 0.8)
dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()Model assumptions
Physiology models
Total Biomass
########## ##########
######### biomass ----
Y<-model.data$biomass
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
summary(full)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season * Light * Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 811
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81107 -0.71811 -0.02691 0.62021 3.10297
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 7.214 2.686
## Residual 59.912 7.740
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.742e+01 2.581e+00 2.116e+01 10.624 6.1e-10
## Seasonwinter -5.409e+00 3.111e+00 1.088e+02 -1.738 0.085
## Light -1.969e-01 2.245e-01 1.107e+02 -0.877 0.382
## DomD 1.722e+00 8.125e+00 1.095e+02 0.212 0.833
## Seasonwinter:Light 7.319e-01 4.225e-01 1.106e+02 1.732 0.086
## Seasonwinter:DomD 2.354e+00 9.431e+00 1.095e+02 0.250 0.803
## Light:DomD 9.707e-04 5.473e-01 1.093e+02 0.002 0.999
## Seasonwinter:Light:DomD -6.400e-02 8.240e-01 1.097e+02 -0.078 0.938
##
## (Intercept) ***
## Seasonwinter .
## Light
## DomD
## Seasonwinter:Light .
## Seasonwinter:DomD
## Light:DomD
## Seasonwinter:Light:DomD
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L Ssn:DD Lgh:DD
## Seasonwintr -0.566
## Light -0.727 0.550
## DomD -0.214 0.190 0.207
## Ssnwntr:Lgh 0.330 -0.818 -0.454 -0.129
## Ssnwntr:DmD 0.182 -0.333 -0.175 -0.865 0.277
## Light:DomD 0.281 -0.236 -0.386 -0.943 0.211 0.815
## Ssnwnt:L:DD -0.172 0.416 0.238 0.636 -0.510 -0.873 -0.671rand(full)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 5.33 1 0.02 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(full, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 184.52 184.52 1 108.7 3.080 0.0821 .
## Light 2.54 2.54 1 105.8 0.042 0.8372
## Dom 43.22 43.22 1 108.3 0.721 0.3975
## Season:Light 232.07 232.07 1 110.5 3.873 0.0516 .
## Season:Dom 3.73 3.73 1 109.5 0.062 0.8034
## Light:Dom 0.27 0.27 1 108.8 0.005 0.9463
## Season:Light:Dom 0.36 0.36 1 109.7 0.006 0.9382
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(full, ~Light:Season)
CLD(posthoc, Letters=letters)## Light Season emmean SE df lower.CL upper.CL .group
## 8.043006 summer 26.70017 2.518114 18.90 21.42787 31.97246 a
## 8.043006 winter 28.09789 1.822194 5.64 23.56837 32.62741 a
##
## Results are averaged over the levels of: Dom
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(full), ylab="biomass", par.strip.text=list(cex=0.7))
Chlorophyll a
######### chltotal--
Y<-model.data$chltot
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
summary(full)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season * Light * Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 419.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.43730 -0.62654 -0.04514 0.46171 3.03177
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.3024 0.5499
## Residual 1.7520 1.3236
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.50685 0.46704 15.02000 11.791 5.40e-09
## Seasonwinter 2.67199 0.53227 108.49000 5.020 2.04e-06
## Light -0.05151 0.03848 110.33000 -1.339 0.1834
## DomD 0.80228 1.39071 109.08000 0.577 0.5652
## Seasonwinter:Light -0.15910 0.07242 110.22000 -2.197 0.0301
## Seasonwinter:DomD -3.36688 1.61430 109.06000 -2.086 0.0393
## Light:DomD -0.10411 0.09366 108.91000 -1.111 0.2688
## Seasonwinter:Light:DomD 0.23738 0.14107 109.24000 1.683 0.0953
##
## (Intercept) ***
## Seasonwinter ***
## Light
## DomD
## Seasonwinter:Light *
## Seasonwinter:DomD *
## Light:DomD
## Seasonwinter:Light:DomD .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L Ssn:DD Lgh:DD
## Seasonwintr -0.533
## Light -0.689 0.547
## DomD -0.201 0.190 0.205
## Ssnwntr:Lgh 0.308 -0.818 -0.447 -0.130
## Ssnwntr:DmD 0.171 -0.333 -0.173 -0.865 0.277
## Light:DomD 0.265 -0.235 -0.385 -0.943 0.210 0.815
## Ssnwnt:L:DD -0.161 0.417 0.235 0.636 -0.510 -0.873 -0.671rand(full)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 7.04 1 0.008 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(full, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 36.79 36.79 1 108.4 20.999 1.24e-05 ***
## Light 18.40 18.40 1 109.0 10.502 0.00158 **
## Dom 10.32 10.32 1 108.1 5.891 0.01687 *
## Season:Light 4.34 4.34 1 110.1 2.475 0.11856
## Season:Dom 7.62 7.62 1 109.1 4.350 0.03934 *
## Light:Dom 0.00 0.00 1 108.5 0.000 0.98600
## Season:Light:Dom 4.96 4.96 1 109.2 2.831 0.09530 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(full, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## summer 5.075014 0.4562797 14.16 4.097437 6.052591 a
## winter 5.738545 0.3463131 4.94 4.845002 6.632088 a
##
## Results are averaged over the levels of: Dom
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(full, ~Dom|Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D 5.057493 0.7525279 62.15 3.553282 6.561704 a
## C 5.092535 0.3386960 4.61 4.199143 5.985928 a
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D 4.992202 0.4220449 10.64 4.059429 5.924976 a
## C 6.484888 0.3917497 7.68 5.574853 7.394923 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(full), ylab="total chlor", par.strip.text=list(cex=0.7))
Chlorophyll per symbiont cell
######### chlcell --
Y<-model.data$chlcell
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 4.8574 21.532 3.5713 -7.1426
## object 10 6.9000 34.691 6.5500 -13.1000 5.9574 4 0.2023summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 13.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7400 -0.5161 -0.0852 0.5185 2.4510
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.02551 0.1597
## Residual 0.05208 0.2282
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.87150 0.09786 5.69000 19.124 2.23e-06 ***
## Seasonwinter 0.07957 0.04766 112.84000 1.669 0.0978 .
## Light -0.03471 0.00525 114.62000 -6.611 1.26e-09 ***
## DomD -0.52057 0.05268 112.94000 -9.881 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.407
## Light -0.478 0.476
## DomD 0.106 -0.261 -0.447rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 23.9 1 1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.1452 0.1452 1 112.84 2.787 0.09779 .
## Light 2.2767 2.2767 1 114.62 43.711 1.261e-09 ***
## Dom 5.0854 5.0854 1 112.94 97.637 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(add), ylab="chla cell", par.strip.text=list(cex=0.7))
#ranef(add)
#rand(add)
#fixef(add)
#sjp.lmer(add, y.offset = .4)Isotope models
host d13C
########## host..d13C --
Y<-model.data$host..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 323.06 339.73 -155.53 311.06
## object 10 326.01 353.80 -153.00 306.01 5.0527 4 0.2819summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 321.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.85036 -0.66017 -0.07034 0.63948 2.03805
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.2418 0.4918
## Residual 0.7654 0.8749
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -16.37006 0.32741 7.43000 -49.998 1.14e-10 ***
## Seasonwinter 0.06996 0.18251 113.27000 0.383 0.702
## Light 0.11577 0.02004 115.00000 5.777 6.62e-08 ***
## DomD -1.12836 0.20172 113.40000 -5.594 1.56e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.465
## Light -0.546 0.474
## DomD 0.119 -0.259 -0.445rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 17.9 1 2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.1124 0.1124 1 113.27 0.147 0.7022
## Light 25.5408 25.5408 1 115.00 33.369 6.615e-08 ***
## Dom 23.9479 23.9479 1 113.40 31.288 1.561e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(add), ylab="d13Chost", par.strip.text=list(cex=0.7))
symbiont d13C
########## symb..d13C --
Y<-model.data$symb..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 340.93 357.61 -164.47 328.93
## object 10 337.87 365.66 -158.94 317.87 11.061 4 0.02589 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(full)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season * Light * Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 339.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.97817 -0.63683 -0.00398 0.73631 2.75134
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.4017 0.6338
## Residual 0.8284 0.9102
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -16.51611 0.41019 7.15000 -40.265 1.09e-09
## Seasonwinter -0.43681 0.36630 108.25000 -1.192 0.2357
## Light 0.13202 0.02657 109.19000 4.968 2.52e-06
## DomD -1.85301 0.95794 108.49000 -1.934 0.0557
## Seasonwinter:Light 0.02116 0.05000 109.11000 0.423 0.6730
## Seasonwinter:DomD 1.66750 1.11193 108.48000 1.500 0.1366
## Light:DomD 0.02498 0.06450 108.41000 0.387 0.6993
## Seasonwinter:Light:DomD -0.03626 0.09720 108.57000 -0.373 0.7099
##
## (Intercept) ***
## Seasonwinter
## Light ***
## DomD .
## Seasonwinter:Light
## Seasonwinter:DomD
## Light:DomD
## Seasonwinter:Light:DomD
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L Ssn:DD Lgh:DD
## Seasonwintr -0.415
## Light -0.542 0.540
## DomD -0.156 0.191 0.202
## Ssnwntr:Lgh 0.235 -0.817 -0.435 -0.131
## Ssnwntr:DmD 0.132 -0.333 -0.170 -0.865 0.278
## Light:DomD 0.207 -0.235 -0.382 -0.943 0.210 0.816
## Ssnwnt:L:DD -0.124 0.417 0.230 0.637 -0.511 -0.874 -0.671rand(full)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 25.5 1 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(full, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.443 0.443 1 108.22 0.535 0.4663
## Light 34.811 34.811 1 110.78 42.021 2.591e-09 ***
## Dom 1.364 1.364 1 108.07 1.646 0.2022
## Season:Light 0.063 0.063 1 109.02 0.076 0.7838
## Season:Dom 1.863 1.863 1 108.48 2.249 0.1366
## Light:Dom 0.029 0.029 1 108.25 0.035 0.8528
## Season:Light:Dom 0.115 0.115 1 108.57 0.139 0.7099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(full), ylab="d13Csymb", par.strip.text=list(cex=0.7))
host-symbiont d13C
######### d13C..host.sym --
Y<-model.data$d13C..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Light +Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 8 28.734 50.967 -6.3670 12.734
## object 10 32.339 60.131 -6.1697 12.339 0.3946 2 0.8209summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 43
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0263 -0.4446 0.0338 0.5208 3.0008
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.01887 0.1374
## Residual 0.06362 0.2522
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.133392 0.096724 8.690000 1.379 0.20233
## Seasonwinter 0.296870 0.092051 110.200000 3.225 0.00166 **
## Light -0.004318 0.006787 112.000000 -0.636 0.52598
## DomD -0.113913 0.087914 110.560000 -1.296 0.19776
## Seasonwinter:Light -0.024891 0.011636 111.280000 -2.139 0.03461 *
## Seasonwinter:DomD -0.288488 0.115782 110.470000 -2.492 0.01420 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L
## Seasonwintr -0.484
## Light -0.587 0.554
## DomD 0.137 -0.111 -0.517
## Ssnwntr:Lgh 0.284 -0.778 -0.484 0.252
## Ssnwntr:DmD -0.086 0.036 0.364 -0.750 -0.357rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 12.4 1 4e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.7091 0.7091 1 110.2 11.147 0.00115 **
## Light 0.2799 0.2799 1 113.0 4.400 0.03817 *
## Dom 1.4788 1.4788 1 111.4 23.246 4.55e-06 ***
## Season:Light 0.2911 0.2911 1 111.3 4.576 0.03461 *
## Season:Dom 0.3950 0.3950 1 110.5 6.208 0.01420 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(full, ~Dom|Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D 0.04352437 0.15115422 40.38 -0.26188100 0.34892974 a
## C 0.09821704 0.07892380 3.97 -0.12150117 0.31793525 a
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D -0.20496599 0.09252625 7.37 -0.42154750 0.01161552 a
## C 0.18935605 0.08747198 5.81 -0.02639938 0.40511148 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(full, ~Dom)
CLD(posthoc, Letters=letters)## Dom emmean SE df lower.CL upper.CL .group
## D -0.08072081 0.1012735 10.39 -0.30523544 0.1437938 a
## C 0.14378654 0.0765190 3.49 -0.08156513 0.3691382 b
##
## Results are averaged over the levels of: Season
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d13C-hs", par.strip.text=list(cex=0.7))
skeleton d13C
####### d13C..skel --
Y<-model.data$d13C..skel
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 311.33 328.01 -149.67 299.33
## object 10 315.61 343.40 -147.81 295.61 3.7238 4 0.4447summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 310.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96903 -0.81931 0.00245 0.67065 2.42406
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.1369 0.3700
## Residual 0.7036 0.8388
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.839026 0.277165 10.130000 -10.243 1.15e-06 ***
## Seasonwinter 0.461515 0.174708 113.770000 2.642 0.00941 **
## Light 0.006456 0.019089 114.340000 0.338 0.73583
## DomD 0.044872 0.193058 113.950000 0.232 0.81662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.523
## Light -0.614 0.471
## DomD 0.132 -0.257 -0.442rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 10.1 1 0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 4.910 4.910 1 113.8 6.978 0.00941 **
## Light 0.080 0.080 1 114.3 0.114 0.73583
## Dom 0.038 0.038 1 114.0 0.054 0.81662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(add, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## summer -2.764663 0.2270645 4.74 -3.358112 -2.171214 a
## winter -2.303148 0.2187764 4.14 -2.902642 -1.703655 b
##
## Results are averaged over the levels of: Dom
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(full), ylab="d13C-skel", par.strip.text=list(cex=0.7))
host d15N
####### host..d15N --
Y<-model.data$host..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Light+ Season:Dom+ Light:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom +
## ..1: (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 9 86.794 111.81 -34.397 68.794
## object 10 88.753 116.55 -34.377 68.753 0.0406 1 0.8404summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom +
## (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 100.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.85472 -0.56530 0.00305 0.63356 2.24611
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.29222 0.5406
## Residual 0.09533 0.3088
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.873694 0.283900 3.570000 17.167 0.000148 ***
## Seasonwinter 0.200690 0.113010 109.040000 1.776 0.078544 .
## Light -0.016754 0.008801 109.290000 -1.904 0.059587 .
## DomD -0.389112 0.250457 109.030000 -1.554 0.123176
## Seasonwinter:Light -0.022145 0.014615 109.190000 -1.515 0.132592
## Seasonwinter:DomD 0.205938 0.183732 109.070000 1.121 0.264808
## Light:DomD 0.032102 0.016224 109.040000 1.979 0.050375 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L Ssn:DD
## Seasonwintr -0.195
## Light -0.259 0.499
## DomD -0.049 -0.107 0.073
## Ssnwntr:Lgh 0.096 -0.773 -0.373 0.293
## Ssnwntr:DmD 0.024 0.069 0.065 -0.822 -0.402
## Light:DomD 0.082 0.067 -0.315 -0.903 -0.210 0.635print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.2754 0.2754 1 109.0 2.888 0.0921 .
## Light 0.5038 0.5038 1 109.7 5.284 0.0234 *
## Dom 0.1139 0.1139 1 109.1 1.195 0.2767
## Season:Light 0.2189 0.2189 1 109.2 2.296 0.1326
## Season:Dom 0.1198 0.1198 1 109.1 1.256 0.2648
## Light:Dom 0.3732 0.3732 1 109.0 3.915 0.0504 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(full), ylab="d15Nhost", par.strip.text=list(cex=0.7))
symbiont d15N
######## symb..d15N --
Y<-model.data$symb..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Light+ Season:Dom+ Light:Dom + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom +
## ..1: (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 9 97.046 122.06 -39.523 79.046
## object 10 98.398 126.19 -39.199 78.398 0.6483 1 0.4207summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + Light:Dom +
## (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 110.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2292 -0.5979 0.1164 0.6161 1.9877
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.4345 0.6592
## Residual 0.1028 0.3206
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.344328 0.341713 3.410000 12.713 0.000546 ***
## Seasonwinter -0.065757 0.117357 109.030000 -0.560 0.576416
## Light -0.004593 0.009141 109.210000 -0.502 0.616330
## DomD -0.390860 0.260091 109.020000 -1.503 0.135786
## Seasonwinter:Light -0.002555 0.015178 109.140000 -0.168 0.866608
## Seasonwinter:DomD 0.488751 0.190804 109.050000 2.562 0.011787 *
## Light:DomD 0.027628 0.016848 109.030000 1.640 0.103925
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L Ssn:DD
## Seasonwintr -0.169
## Light -0.224 0.499
## DomD -0.042 -0.107 0.073
## Ssnwntr:Lgh 0.083 -0.773 -0.372 0.293
## Ssnwntr:DmD 0.021 0.069 0.065 -0.822 -0.402
## Light:DomD 0.071 0.067 -0.315 -0.903 -0.210 0.635print(anova(add, type=2), digits=3)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.057 0.057 1 109 0.56 0.457
## Light 0.001 0.001 1 110 0.01 0.937
## Dom 0.124 0.124 1 109 1.21 0.274
## Season:Light 0.003 0.003 1 109 0.03 0.867
## Season:Dom 0.675 0.675 1 109 6.56 0.012 *
## Light:Dom 0.276 0.276 1 109 2.69 0.104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(full, ~Season*Dom)
CLD(posthoc, Letters=letters)## Season Dom emmean SE df lower.CL upper.CL .group
## winter C 4.206022 0.3338187 3.20 3.180359 5.231685 a
## summer D 4.210158 0.3683973 4.74 3.247190 5.173125 ab
## summer C 4.306919 0.3304207 3.07 3.269373 5.344465 a
## winter D 4.546253 0.3360060 3.29 3.527659 5.564846 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05plot(allEffects(full), ylab="d15Nsymb", par.strip.text=list(cex=0.7))
host-symbiont d15N
####### d15N..host.sym --
Y<-model.data$d15N..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 7 1.6903 21.144 6.1548 -12.310
## object 10 3.5687 31.360 8.2156 -16.431 4.1216 3 0.2486summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 12
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.80715 -0.73247 0.06926 0.71888 2.22092
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.01318 0.1148
## Residual 0.05131 0.2265
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.555024 0.080132 8.180000 6.926 0.000109 ***
## Seasonwinter 0.155035 0.051874 112.000000 2.989 0.003444 **
## Light -0.015260 0.005323 113.890000 -2.867 0.004942 **
## DomD 0.090965 0.076384 112.290000 1.191 0.236207
## Seasonwinter:DomD -0.380440 0.097092 111.940000 -3.918 0.000154 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD
## Seasonwintr -0.454
## Light -0.556 0.322
## DomD 0.072 0.141 -0.465
## Ssnwntr:DmD 0.018 -0.414 0.233 -0.730print(anova(add, type=2), digits=3)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.107 0.107 1 112 2.08 0.15217
## Light 0.422 0.422 1 114 8.22 0.00494 **
## Dom 0.308 0.308 1 113 6.01 0.01580 *
## Season:Dom 0.788 0.788 1 112 15.35 0.00015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(full, ~Season*Dom)
CLD(posthoc, Letters=letters)## Season Dom emmean SE df lower.CL upper.CL .group
## winter D 0.2892355 0.08014858 7.79 0.1035443 0.4749266 a
## summer D 0.4009953 0.13294032 43.72 0.1330241 0.6689665 ab
## summer C 0.4314180 0.06772961 4.06 0.2444074 0.6184287 ab
## winter C 0.5608237 0.07554987 6.06 0.3763859 0.7452615 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d15N-hs", par.strip.text=list(cex=0.7))
host C:N
###### host..C.N --
Y<-model.data$host..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 -271.49 -254.82 141.75 -283.49
## object 10 -266.00 -238.21 143.00 -286.00 2.5084 4 0.6431summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: -252.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8291 -0.7150 -0.1038 0.6546 3.2395
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0008875 0.02979
## Residual 0.0052739 0.07262
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.946e+00 2.324e-02 1.054e+01 83.720 4.44e-16 ***
## Seasonwinter 4.134e-03 1.512e-02 1.139e+02 0.274 0.785
## Light 1.116e-03 1.649e-03 1.137e+02 0.677 0.500
## DomD 1.002e-02 1.670e-02 1.141e+02 0.600 0.550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.539
## Light -0.633 0.470
## DomD 0.135 -0.256 -0.441print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.0003945 0.0003945 1 113.9 0.0748 0.785
## Light 0.0024172 0.0024172 1 113.7 0.4583 0.500
## Dom 0.0018996 0.0018996 1 114.1 0.3602 0.550plot(allEffects(add), ylab="C:N host", par.strip.text=list(cex=0.7))
symbiont C:N
####### symb..C.N --
Y<-model.data$symb..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use full model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 -210.60 -193.98 111.30 -222.60
## object 10 -205.91 -178.21 112.96 -225.91 3.3076 4 0.5077summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: -192.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8681 -0.6308 -0.1693 0.6705 2.9795
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 4.491e-18 2.119e-09
## Residual 9.188e-03 9.585e-02
## Number of obs: 118, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.053e+00 2.219e-02 1.140e+02 92.531 <2e-16 ***
## Seasonwinter -2.958e-02 1.956e-02 1.140e+02 -1.513 0.133
## Light 4.094e-06 1.982e-03 1.140e+02 0.002 0.998
## DomD 5.402e-03 2.150e-02 1.140e+02 0.251 0.802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.680
## Light -0.800 0.426
## DomD 0.116 -0.226 -0.400print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.0210230 0.0210230 1 114 2.28813 0.1331
## Light 0.0000000 0.0000000 1 114 0.00000 0.9984
## Dom 0.0005801 0.0005801 1 114 0.06314 0.8021plot(allEffects(full), ylab="C:N symb", par.strip.text=list(cex=0.7))
Figures
Physiology
#########################
# Physiology
df.fig<-data.trim
# figure layout
layout(matrix(c(1,1,2,3), 2,2, byrow=TRUE))
par(mar=c(5,4.5,2,2))
##################
##################
## season relationship
data.summer<-df.fig[df.fig$Season=="summer", ]
data.winter<-df.fig[df.fig$Season=="winter", ]
## Site relationship
F8.R10.df<-df.fig[df.fig$Reef=="F8-R10", ]
F1.R46.df<-df.fig[df.fig$Reef=="F1-R46", ]
HIMB.df<-df.fig[df.fig$Reef=="HIMB", ]
R42.df<-df.fig[df.fig$Reef=="R42", ]
## Symbiont relationship
C.sum.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="summer") ,]
C.win.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="winter") ,]; C.win.df<-na.omit(C.win.df)
D.sum.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="summer") ,]; D.sum.df<-na.omit(D.sum.df)
D.win.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="winter") ,]; D.win.df<-na.omit(D.win.df)
##################
# Fig: chlorophyll (total) over season and depth
##################
plot(chltot~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")),
main="chl: Depth and Season",
xlim=c(0,11), ylim=c(1,17))
abline((lm(chltot~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points(chltot~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm(chltot~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(data.winter$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="dodgerblue3", main="chl: Seasons",
xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")))
lines(density(data.summer$chltot), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(F1.R46.df$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="tomato2", main="chl: Sites",
xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), ylab="Density")
lines(density(R42.df$chltot), col="skyblue3")
lines(density(F8.R10.df$chltot), col="springgreen4")
lines(density(HIMB.df$chltot), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.chl.pdf", width=9, height=7)
##################
# Fig: chlorophyll (total) over season and depth
##################
plot(chltot~Light, data=C.sum.df,col="tomato2", pch=16, cex=0.7,
xlab="Light (DLI)", ylab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")),
main="chl: Light and Season",
xlim=c(0,20), ylim=c(1,12))
abline((lm(chltot~Light, data=C.sum.df)), col='tomato2', lwd=2)
points(chltot~Light, data=C.win.df, col="dodgerblue3", pch=16, cex=0.7)
abline((lm(chltot~Light, data=C.win.df)), col='dodgerblue3', lwd=2)
points(chltot~Light, data=D.sum.df, col="mediumseagreen", pch=16, cex=0.7)
abline((lm(chltot~Light, data=D.sum.df)), col='mediumseagreen', lwd=2)
points(chltot~Light, data=D.win.df, col="orchid", pch=16, cex=0.7)
abline((lm(chltot~Light, data=D.win.df)), col='orange', lwd=2)
legend("topright", c("C-sum", "C-win", "D-sum", "D-win"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3","mediumseagreen", "orange"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.05))
plot(density(C.sum.df$chltot), col="tomato2", pch=16, cex=0.7,
xlab="Light (DLI)", ylab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")),
main="chl: Light and Season", ylim=c(0, 0.4), xlim=c(0,15), lwd=2)
lines(density(C.win.df$chltot), col="dodgerblue3", lwd=2)
lines(density(D.sum.df$chltot), col="orange", lwd=2)
lines(density(D.win.df$chltot), col="mediumseagreen", lwd=2)
legend("topright", c("C-sum", "C-win", "D-sum", "D-win"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3", "orange", "mediumseagreen"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.05))
###### density plot: by season
plot(density(data.winter$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="dodgerblue3", main="chl: Seasons",
xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")))
lines(density(data.summer$chltot), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(F1.R46.df$chltot), ylim=c(0,0.3), xlim=c(0, 18), col="tomato2", main="chl: Sites",
xlab=expression(paste("chlorophyll", ~(mu*g~cm^-2), sep="")), ylab="Density")
lines(density(R42.df$chltot), col="skyblue3")
lines(density(F8.R10.df$chltot), col="springgreen4")
lines(density(HIMB.df$chltot), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.chl.pdf", width=9, height=7)
##################
# Fig: chlorophyll/cell over season and depth
##################
plot(chlcell~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)",
ylab=expression(paste("pg chlorophyll cell"^-1, sep="")), main="chl/cell: Depth and Season",
xlim=c(0,10), ylim=c(0,15))
abline((lm(chlcell~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points(chlcell~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm(chlcell~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(data.winter$chlcell), ylim=c(0,0.4), xlim=c(0, 15), col="dodgerblue3", main="chl/cell: Seasons", xlab=expression(paste("pg chlorophyll cell"^-1, sep="")))
lines(density(data.summer$chlcell), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(F1.R46.df$chlcell), ylim=c(0,0.4), xlim=c(0,15), col="tomato2", main="chl/cell: Sites", xlab=expression(paste("pg chlorophyll cell"^-1, sep="")))
lines(density(R42.df$chlcell), col="skyblue3")
lines(density(F8.R10.df$chlcell), col="springgreen4")
lines(density(HIMB.df$chlcell), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.chlcell.pdf", width=9, height=7)
##################
# Fig: symbionts over season and depth
##################
plot((zoox/10^6)~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab=(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))), xlim=c(0,10), ylim=c(0, 6),
main= "zoox: Depth and Season")
abline((lm((zoox/10^6)~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points((zoox/10^6)~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm((zoox/10^6)~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(data.summer$zoox/10^6), col="tomato2", main="zoox: Seasons",
xlab=(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))), xlim=c(0, 7))
lines(density(data.winter$zoox/10^6), col="dodgerblue3")
legend('topright', c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(F1.R46.df$zoox/10^6), col="tomato2", main="zoox: Sites",
xlab=(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))), xlim=c(0, 7))
lines(density(R42.df$zoox/10^6), col="skyblue3")
lines(density(F8.R10.df$zoox/10^6), col="springgreen4")
lines(density(HIMB.df$zoox/10^6), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.zoox.pdf", width=9, height=7)
##################
# Fig: biomass over season and depth
##################
plot(biomass~newDepth, data=data.winter,col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste("Total biomass", ~(mg~cm^-2), sep="")),
xlim=c(0,10), ylim=c(5,60),
main="biomass: Depth and Season")
abline((lm(biomass~newDepth, data=data.winter)), col='dodgerblue3', lwd=2)
points(biomass~newDepth, data=data.summer, col="tomato2", pch=16, cex=0.7)
abline((lm(biomass~newDepth, data=data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(data.winter$biomass), ylim=c(0,0.08), xlim=c(0, 70), col="dodgerblue3", main="biomass: Seasons", xlab=expression(paste("Total biomass", ~(mg~cm^-2), sep="")))
lines(density(data.summer$biomass), col="tomato2")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), cex=1, y.intersp = 0.3, x.intersp = 0.4, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(F1.R46.df$biomass), ylim=c(0,0.1), xlim=c(0,70), col="tomato2", main="biomass: Sites",
xlab=expression(paste("Total biomass", ~(mg~cm^-2), sep="")))
lines(density(R42.df$biomass), col="skyblue3")
lines(density(F8.R10.df$biomass), col="springgreen4")
lines(density(HIMB.df$biomass), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2),
y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/physiology/PanKB.biomass.pdf", width=9, height=7)Isotopes
### ISOTOPES ##
###############
##################
# Fig: d13C host over season and depth
##################
plot(host..d13C~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)")),
xlim=c(0,10), ylim=c(-20,-10),
main= "d13C host: Depth and Season")
abline(lm(host..d13C~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(host..d13C~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(host..d13C~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$host..d13C)), col="tomato2", main="d13C-host: Seasons",
xlab=expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)")), xlim=c(-20, -8))
lines(density(na.omit(data.winter$host..d13C)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$host..d13C)), col="tomato2", main="d13C-host: Sites",
xlab=expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)")), ylim=c(0,0.4), xlim=c(-20, -8))
lines(density(na.omit(R42.df$host..d13C)), col="skyblue3")
lines(density(na.omit(F8.R10.df$host..d13C)), col="springgreen4")
lines(density(na.omit(HIMB.df$host..d13C)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d13C-host.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: d13C symb over season and depth
##################
plot(symb..d13C~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)")),
xlim=c(0,10), ylim=c(-20,-10),
main= "d13C symb: Depth and Season")
abline(lm(symb..d13C~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(symb..d13C~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(symb..d13C~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$symb..d13C)), col="tomato2", main="d13C-symb: Seasons",
xlab=expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)")), xlim=c(-20, -8))
lines(density(na.omit(data.winter$symb..d13C)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$symb..d13C)), col="tomato2", main="d13C-symb: Sites",
xlab=expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)")), ylim=c(0,0.4), xlim=c(-20, -8))
lines(density(na.omit(R42.df$symb..d13C)), col="skyblue3")
lines(density(na.omit(F8.R10.df$symb..d13C)), col="springgreen4")
lines(density(na.omit(HIMB.df$symb..d13C)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d13C-symb.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: d13C skeleton over season and depth
##################
plot(d13C..skel~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[Skel], " (\u2030, V-PDB)")),
xlim=c(0,10), ylim=c(-6,2),
main= "d13C skel: Depth and Season")
abline(lm(d13C..skel~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(d13C..skel~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(d13C..skel~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$d13C..skel)), col="tomato2", main="d13C-skel: Seasons",
xlab=expression(paste(delta^{13}, C[Skel], " (\u2030, V-PDB)")), ylim=c(0, 0.5), xlim=c(-8, 6))
lines(density(na.omit(data.winter$d13C..skel)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$d13C..skel)), col="tomato2", main="d13C-skel: Sites",
xlab=expression(paste(delta^{13}, C[Skel], " (\u2030, V-PDB)")), ylim=c(0,0.5), xlim=c(-8, 6))
lines(density(na.omit(R42.df$d13C..skel)), col="skyblue3")
lines(density(na.omit(F8.R10.df$d13C..skel)), col="springgreen4")
lines(density(na.omit(HIMB.df$d13C..skel)), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d13C-skel.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: d13C host-symb over season and depth
##################
plot(d13C..host.sym~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)")),
xlim=c(0,10), ylim=c(-2.5,3),
main= "d13C h-s: Depth and Season")
abline(lm(d13C..host.sym~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(d13C..host.sym~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(d13C..host.sym~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$d13C..host.sym)), col="tomato2", main="d13C h-s: Seasons",
xlab=expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)")), ylim=c(0,1.5), xlim=c(-3, 5))
lines(density(na.omit(data.winter$d13C..host.sym)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$d13C..host.sym)), col="tomato2", main="d13C h-s: Sites",
expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)")), ylim=c(0,1.7), xlim=c(-3,5))
lines(density(na.omit(R42.df$d13C..host.sym)), col="skyblue3")
lines(density(na.omit(F8.R10.df$d13C..host.sym)), col="springgreen4")
lines(density(na.omit(HIMB.df$d13C..host.sym)), col="purple")
legend("topright",c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d13Ch-s.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: d15N host over season and depth
##################
plot(host..d15N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste(delta^{15}, N[H], " (\u2030, air)")), xlim=c(0,10), ylim=c(2,8),
main= "d15N host: Depth and Season")
abline(lm(host..d15N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(host..d15N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(host..d15N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$host..d15N)), col="tomato2", main="d15N-host: Seasons",
xlab=expression(paste(delta^{15}, N[H], " (\u2030, air)")), xlim=c(0, 12), ylim=c(0, 1))
lines(density(na.omit(data.winter$host..d15N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$host..d15N)), col="tomato2", main="d15N-host: Sites", xlab=expression(paste(delta^{15}, N[H], " (\u2030, air)")), ylim=c(0,2), xlim=c(0, 12))
lines(density(na.omit(R42.df$host..d15N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$host..d15N)), col="springgreen4")
lines(density(na.omit(HIMB.df$host..d15N)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d15N-host.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: d15N symb over season and depth
##################
plot(symb..d15N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab =expression(paste(delta^{15}, N[S], " (\u2030, air)")),
xlim=c(0,10), ylim=c(2,7),
main= "d15N symb: Depth and Season")
abline(lm(symb..d15N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(symb..d15N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(symb..d15N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$symb..d15N)), col="tomato2", main="d15N-symb: Seasons",
xlab= expression(paste(delta^{15}, N[S], " (\u2030, air)")), xlim=c(0, 9), ylim=c(0, 0.9))
lines(density(na.omit(data.winter$symb..d15N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$symb..d15N)), col="tomato2", main="d15N-symb: Sites",
xlab= expression(paste(delta^{15}, N[S], " (\u2030, air)")), xlim=c(0, 9), ylim=c(0,1.8))
lines(density(na.omit(R42.df$symb..d15N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$symb..d15N)), col="springgreen4")
lines(density(na.omit(HIMB.df$symb..d15N)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d15N-symb.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: d15N host.symb over season and depth
##################
plot(d15N..host.sym~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste(delta^{15}, N[H-S], " (\u2030, air)")),
xlim=c(0,10), ylim=c(-1,2),
main= "d15N h-s: Depth and Season")
abline(lm(d15N..host.sym~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(d15N..host.sym~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(d15N..host.sym~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$d15N..host.sym)), col="tomato2", main="d15N h-s: Seasons",
xlab=expression(paste(delta^{15}, N[H-S], " (\u2030, air)")), xlim=c(-1, 3), ylim=c(0,2))
lines(density(na.omit(data.winter$d15N..host.sym)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.3, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$d15N..host.sym)), col="tomato2", main="d15N h-s: Sites", xlab=expression(paste(delta^{15}, N[H-S], " (\u2030, air)")), ylim=c(0,2), xlim=c(-1,3))
lines(density(na.omit(R42.df$d15N..host.sym)), col="skyblue3")
lines(density(na.omit(F8.R10.df$d15N..host.sym)), col="springgreen4")
lines(density(na.omit(HIMB.df$d15N..host.sym)), col="purple")
legend("topright" ,c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.d15Nh-s.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: C.N host season and depth
##################
plot(host..C.N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste("C:N"[H])),
xlim=c(0,10), ylim=c(5,10),
main= "C:N-host Depth and Season")
abline(lm(host..C.N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(host..C.N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(host..C.N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$host..C.N)), col="tomato2", main="C:N-host Depth and Season",
xlab=expression(paste("C:N"[H])), ylim=c(0, 1), xlim=c(4, 10))
lines(density(na.omit(data.winter$host..C.N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$host..C.N)), col="tomato2", main="C:N-host Depth and Season",
xlab=expression(paste("C:N"[H])), ylim=c(0,1.5), xlim=c(4, 10))
lines(density(na.omit(R42.df$host..C.N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$host..C.N)), col="springgreen4")
lines(density(na.omit(HIMB.df$host..C.N)), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.C.Nhost.pdf", width=9, height=7, encod="MacRoman")
##################
# Fig: C.N symb season and depth
##################
plot(symb..C.N~newDepth, data=na.omit(data.winter),col="dodgerblue3", pch=16, cex=0.7,
xlab="Depth (m)", ylab = expression(paste("C:N"[S])),
xlim=c(0,10), ylim=c(5,14),
main= "C:N-symb Depth and Season")
abline(lm(symb..C.N~newDepth, data=na.omit(data.winter)), col='dodgerblue3', lwd=2)
points(symb..C.N~newDepth, data=na.omit(data.summer), col="tomato2", pch=16, cex=0.7)
abline(lm(symb..C.N~newDepth, data=na.omit(data.summer)), col='tomato2', lwd=2)
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.08, -0.1))
###### density plot: by season
plot(density(na.omit(data.summer$symb..C.N)), col="tomato2", main="C:N-symb Depth and Season",
xlab=expression(paste("C:N"[S])), ylim=c(0, 0.8), xlim=c(5, 15))
lines(density(na.omit(data.winter$symb..C.N)), col="dodgerblue3")
legend("topright", c("summer", "winter"), lty=c(1,1), lwd=c(2,2),
col=c("tomato2", "dodgerblue3"), y.intersp = 0.3, x.intersp = 0.4, cex=1, bty="n", inset=c(-0.35, -0.1))
###### density plot: by Site
plot(density(na.omit(F1.R46.df$symb..C.N)), col="tomato2", main="C:N-symb Depth and Season",
xlab=expression(paste("C:N"[S])), ylim=c(0,0.8), xlim=c(5, 15))
lines(density(na.omit(R42.df$symb..C.N)), col="skyblue3")
lines(density(na.omit(F8.R10.df$symb..C.N)), col="springgreen4")
lines(density(na.omit(HIMB.df$symb..C.N)), col="purple")
legend("topright", c("Fringe-Reef 46", "Patch-Reef 42", "Fringe-Reef 10", "HIMB"), lty=c(1,1), lwd=c(2,2), y.intersp = 0.3, x.intersp = 0.4,
col=c("sienna1", "skyblue3", "springgreen4", "purple"), cex=0.8, bty='n', inset=c(-0.5, -0.1))
dev.copy(pdf, "figures/isotope/PanKB.C.Nsym.pdf", width=9, height=7, encod="MacRoman")